Available  online  at  www.sciencedirect.com 


ELSEVIER 


%V  ScienceDirect 


Journal  of  Power  Sources  160  (2006)  111  1-1 121 


JOURNAL  OF 


www.elsevier.com /locate /jpowsour 


Numerical  simulation  of  a  mini  PEMFC  stack 

Zhixiang  Liua,  Zongqiang  Maoa’*,  Cheng  Wanga,  Weilin  Zhuge,  Yangjun  Zhangb 

a  Institute  of  Nuclear  and  New  Energy  Technology,  Tsinghua  University,  Beijing  100084,  PR  China 
b  Department  of  Automotive  Engineering  and  the  State  Key  Laboratory  of  Automotive  Safety  and  Energy,  Tsinghua  University,  Beijing  100084,  PR  China 

Received  9  January  2006;  received  in  revised  form  28  February  2006;  accepted  1  March  2006 

Available  online  17  April  2006 


Abstract 

Fuel  cell  modeling  and  simulation  has  aroused  much  attention  recently  because  it  can  probe  transport  and  reaction  mechanism.  In  this  paper,  a 
computational  fuel  cell  dynamics  (CFCD)  method  was  applied  to  simulate  a  proton  exchange  membrane  fuel  cell  (PEMFC)  stack  for  the  first  time. 
The  air  cooling  mini  fuel  cell  stack  consisted  of  six  cells,  in  which  the  active  area  was  8  cm2  (2  cm  x  4  cm).  With  reasonable  simplification,  the 
computational  elements  were  effectively  reduced  and  allowed  a  simulation  which  could  be  conducted  on  a  personal  computer  without  large-scale 
parallel  computation.  The  results  indicated  that  the  temperature  gradient  inside  the  fuel  cell  stack  was  determined  by  the  flow  rate  of  the  cooling 
air.  If  the  air  flow  rate  is  too  low,  the  stack  could  not  be  effectively  cooled  and  the  temperature  will  rise  to  a  range  that  might  cause  unstable  stack 
operation. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  fuel  cells  (PEMFCs)  are  promis¬ 
ing  power  sources  for  stationary  power  stations,  electric  vehi¬ 
cles  and  portable  devices  [1].  Almost  all  key  players  in  auto¬ 
mobile  manufacturing  are  now  engaged  in  PEMFC  vehicle 
development. 

Since  the  pioneering  work  of  Springer  et  al.  [2]  and  Bernadi 
and  Verbrugge  [3],  fuel  cell  modeling  has  attracted  the  research 
interest  of  many  groups,  because  modeling  could  provide 
detailed  transport  and  reaction  information  inside  the  fuel  cells 
and  guide  the  designs  of  fuel  cell  engineers.  Costamagna  and 
Srinivasan  [4],  Yao  et  al.  [5],  Wang  [6],  Biyikoglu  [7]  and  Faghri 
and  Guo  [8]  have  given  good  reviews  on  fuel  cell  modeling. 
Since  Gurau  et  al.  [9]  introduced  computational  fluid  dynam¬ 
ics  (CFD)  method  into  fuel  cell  modeling,  CFD  based  fuel  cell 
modeling  has  achieved  great  success  and  commercial  CFD  pack¬ 
ages  are  available  now.  Wang  [6]  termed  CFD  coupled  with 
fuel  cell  transport  and  reaction  modeling  as  computational  fuel 
cell  dynamics  (CFCD).  To  date,  many  papers  on  CFCD  have 
been  published  for  single  cell  modeling  and  simulation,  some 
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of  which  focused  on  the  transport  and  reaction  in  a  segment  of 
the  flow  channel  [10-12],  and  some  others  produced  flow  pat¬ 
terns  [13,14].  There  was  no  CFCD  modeling  or  simulation  work 
available  in  public  documents  for  the  fuel  cell  stack.  One  of  the 
main  difficulties  in  stack  CFCD  is  the  amount  of  computation 
needed.  Including  the  flow  channels  in  each  cell  plate  of  the 
stack,  the  computational  gridpoints  might  be  in  the  billions.  In 
recent  years,  large-scale  parallel  computation  has  been  success¬ 
fully  applied  [15-18]  in  big  area  single  fuel  cell  simulations. 
Wang  and  Wang  [18]  reported  an  ultra  large-scale  simulation 
with  23.5  million  gridpoints,  which  was  the  largest  scale  com¬ 
putation  reported.  However,  CFCD  for  fuel  cell  stack  simulation 
is  still  difficult  without  simplification.  In  this  paper,  we  present 
CFCD  simulation  results  for  a  small  PEMFC  stack  with  reason¬ 
able  simplification,  which  is  the  first  report  of  a  3D  CFD  based 
fuel  cell  stack  simulation. 

2.  Experimental 

Fig.  1  shows  a  photograph  of  the  experimental  mini  fuel  cell 
stack  with  nine  cells.  E-Teck  Pt/C  catalyst  (0.4  mg  cm-2),  Toray 
carbon  paper  (190  [xm  thick)  and  homemade  self-humidifying 
membrane  (50  |xm  thick)  were  used  to  prepare  the  membrane 
electrolyte  assemblies  (MEAs).  The  active  area  was  8  cm2 
(2  cm  x  4  cm).  Graphite  bipolar  plates  with  straight  air  channels 
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the  hydrogen  and  air  manifolds  were  also  included  in  the  com¬ 
putation  domain. 

3.2.  Model  equations 

Governing  equations  of  the  model  employed  in  this  work 
include  conservation  of  mass,  momentum,  energy,  species  and 
current  [19,20]: 


Fig.  1.  Experimental  mini  fuel  cell  stack. 


Mass  conservation  equation: 

2-  Op)  +  v  •  (spU)  =  o  (l) 

ot 

Momentum  conservation  equation: 


(width:  2  mm,  depth:  2  mm)  and  a  serpentine  hydrogen  channel 
(width:  2  mm,  depth:  0.8  mm)  were  assembled  in  the  stack.  A 
micro-air  fan  was  assembled  in  the  stack  to  deliver  air  to  the 
cells. 

During  the  performance  measurement,  the  stack  was  placed 
in  a  humidistat  and  an  Arbin  fuel  cell  test  stand  was  used  to 
collect  performance  data.  The  temperature  and  relative  humidity 
of  the  air  inside  humidistat  were  adjusted  to  be  300  K  and  80% 
relatively. 

3.  Model  description 

3.1.  Model  domain 

The  model  domain  is  shown  in  Fig.  2.  For  the  purpose  of 
limiting  the  number  of  gridpoints  to  make  the  computation  pos¬ 
sible  on  a  personal  computer,  six  cells  were  simulated  instead 
of  nine  cells  as  in  the  experimental  stack.  As  shown  in  Fig.  2, 
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Fig.  2.  Model  domain  of  the  fuel  cell  stack. 


Fig.  3.  (a  and  b)  Simplification  of  the  gas  flow  fields. 
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Energy  conservation  equation: 


3.3.  Simplification  and  grids 


d  dp 

—  (sph)  +  V  •  ( spUh )  =  V  •  a  +  s  — 
dt  H  d  t 


Species  conservation  equation: 
d 


l  X  l 

ho  H - b 

o 


(3) 


dt 


(epY a)  +  V  .  (spUYi)  =  V  •  /*■  +  (hi 


(4) 


Current  conservation: 


V  •  (<TfV0f)  =  •  (^sV^s)  =  JT 


(5) 


And  the  B utter- Volmer  equation  for  the  electrochemical  reac¬ 
tion: 


J  T  =  JO 


a^F 


n  N 


exp 


RT 


0 


exp 


(6) 


More  detail  introduction  of  the  model  could  be  found  in  [20]. 
Liquid  water  is  not  taken  into  account  and  ideal  gas  law  is 
applied. 


For  the  purpose  of  reducing  grids  number,  the  flow  fields 
were  simplified  to  be  porous  media.  As  shown  in  Fig.  3,  the 
bipolar  plate  with  air  and  hydrogen  flow  fields  was  simplified 
to  be  two  porous  layers  separated  by  a  non-permeable  plate. 
The  porosity  of  the  porous  anode  and  cathode  flow  fields  was 
set  to  be  0.67  because  the  flow  channels  took  2/3  of  the  total 
volume.  Accordingly,  the  volume  related  parameters  such  as 
electric  conductivity  and  thermal  conductivity  were  set  to  be 
1/3  of  the  graphite  plate.  Tortuosity  of  the  cathode  flow  field 
was  set  to  be  1.0  because  the  cathode  flow  channels  were 
straight,  and  1.5  was  set  for  the  anode.  To  simulate  the  pres¬ 
sure  drop  in  the  porous  flow  field,  a  section  of  fuel  cell  MEA 
with  10  cm  straight  channels  (as  shown  in  Fig.  3(b))  was  mod¬ 
eled  to  collect  the  pressure  drop  data  with  different  current 
density  discharges.  With  the  pressure  drop  data,  Darcy’s  law 
was  applied  to  calculate  the  permeability  of  the  porous  flow 
field: 


k  = 


fjL  Qy 

A p  ’  ~S~ 


(7) 


Y 

X 


Fig.  4.  Grids  of  the  computational  domain. 
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where  k  is  the  permeability  of  the  equivalent  porous  media  (m2); 
li  the  viscosity  of  air  (Pas);  L  the  channel  length  (m);  A p  the 
pressure  drop  in  the  flow  channel  (Pa);  Qy  the  volume  flux  of 
air  (m3  s_1);  and  S  is  the  section  area  of  the  flow  channel  (m2). 
The  calculated  data  did  not  vary  quite  much  and  an  average  value 
2.4  x  10-8  m2  was  applied  for  the  permeability  of  the  anode  and 
cathode  porous  flow  fields. 

The  computation  domain  was  discreted  to  generate  struc¬ 
tural  grids  with  CFD-Micromesh®.  The  grids  of  the  domain 
are  shown  in  Fig.  4  in  both  xy  and  yz  profiles.  In  xy  profile, 
the  plane  was  evenly  discreted  to  cells  with  1  mm2  area.  In 
yz  profile,  the  grids  were  locally  exponential  distributed  and 
there  were  150  elements  in  z  direction.  With  the  flow  channel 
been  simplified  to  porous  media,  the  grid  number  in  xy  pro¬ 
file  could  be  reduced  without  taking  the  channel  and  rib  width 
into  account.  The  total  grids  were  less  than  200,000.  With¬ 
out  the  simplification,  the  grid  number  should  be  at  least  four 
times  greater,  which  is  beyond  the  computation  capacity  of  a 
microcomputer. 

3.4.  Solution  method 

Parameters  in  the  model  are  listed  in  Table  1.  Anode  inlet 
velocity  is  fixed  to  be  0.5  ms-1,  and  the  cathode  inlet  velocity 
varies  with  different  operation  conditions.  Compared  with  the 
convective  heat  exchanging  in  the  flow  fields,  heat  dissipation 
by  natural  convection  and  radiation  from  the  stack  surface  is 
relative  small.  Additionally,  the  end  plate  of  the  experimental 
stack  is  plastic  and  thermal  conductivity  is  very  small.  So  the 
surface  walls  of  the  computational  domain  were  assumed  to  be 
adiabatic. 

The  problem  was  solved  with  CFD-ACE+  solver.  A  personal 
computer  with  Intel  Pentium  IV  2.4  GHz  processor  and  1.0  GB 


Fig.  5.  I-V  and  I-P  curves  of  the  stack. 
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Fig.  6.  Potential  distribution  across  the  fuel  cell  stack  at  4  A  discharge. 


Table  1 

Parameters  in  the  model 


Parameters 

Value 

Parameters 

Value 

Thickness 

Cathode  flow  field 

2mm 

Thermal  conductivity 

Plate 

129  W  (m  K)-1 

Anode  flow  field 

0.5  mm 

Flow  field 

43  W  (m  K)-1 

Plate 

0.5  mm 

GDL  and  CL 

21  W(mK)-1 

Gas  diffusion  layer  (GDL) 

200  (xm 

PEM 

0.67  W(mK)-1 

Catalyst  layer  (CL) 

25  [xm 

Electric  conductivity 

Proton  conducting  membrane  (PEM) 

50  pan 

Plate 

30008  m-1 

Porosity 

Flow  field 

1000S  m-1 

Flow  field 

0.67 

GDL  and  CL 

12508  m-1 

GDL 

0.6 

PEM 

1.8  x  10-20  S  m-1 

CL 

0.5 

Ionic  conductivity 

PEM 

0.28 

CL 

4.2  8  m-1 

Tortuosity 

PEM 

Eq.  from  [2] 

GDL  and  CL 

1.5 

Exchange  current  density 

PEM 

10 

H2  oxidation 

1.0  x  109Am“3 

Cathode  flow  field 

1.0 

O2  reduction 

2.0  x  105  Am-3 

Anode  flow  field 

1.5 

Operation  condition 

Permeability 

F^/air  temperature 

300  K 

Flow  field 

2.4  x  10-8  m2 

F^/air  relative  humidity 

80% 

GDL  and  CL 

1.76  x  10“ 11  m2 

Outlet  F^/air  pressure 

0.1  MPa 

PEM 

1.8  x  10“18m2 
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RAM  was  adopted  for  the  computation.  The  calculation  con¬ 
verged  after  around  2500  iterations  and  took  around  30  h  for 
each  case. 

4.  Results  and  discussions 

Fig.  5  compares  the  averaged  cell  I-V  and  I-P  curves  of  the 
stack  from  simulation  and  experiments.  In  many  models,  the 
thermodynamic  potential  was  applied  to  calculate  the  open  cir¬ 
cuit  voltage  (OCV)  of  the  fuel  cell  [21,22];  some  authors  used 
empirical  equations  to  calculate  this  parameter  [23].  OCV  cal¬ 
culated  from  these  methods  is  higher  than  1 . 1  V  for  normal  fuel 
cell  operation  conditions.  However,  normally  cell  OCV  is  not 
higher  than  1.0  V  in  experiments.  In  this  model,  0.92  V  aver¬ 
age  cell  OCV  from  experiments  is  used.  In  Fig.  5,  it  can  be 
seen  that  the  simulation  curves  are  in  good  accordance  with  the 
experimental  curves  when  the  current  density  is  relatively  low. 
In  the  high  current  density  region,  the  experimental  I-V  curve 


deviates  greatly  from  linear  relation,  while  the  simulation  curve 
shows  good  linear  behavior.  The  performance  curves  show  a 
great  difference  in  the  high  current  density  region.  One  of  the 
reasons  for  this  difference  is  that  current  density  beneath  the  ribs 


z  x 

(b) 


Fig.  7.  Velocity  and  pressure  distribution  in  the  stack  at  4  A  discharge,  (a)  Air  Fig.  8.  Temperature  distribution  of  the  fuel  cell  stack  in  xz  and  yx  section  planes 
velocity  (ms-1)  and  (b)  total  pressure  (Pa).  (a)  and  in  xy  section  plane  (b)  with  5ms-1  air  inlet  velocity  and  4  A  discharge. 
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is  much  lower  than  that  beneath  the  channel  when  the  fuel  cell 
operates  in  the  low  voltage  region  [22].  But  in  the  simulation, 
because  the  flow  field  is  simplified  to  be  a  porous  media,  the 
influence  of  the  rib  is  neglected.  Another  reason  might  be  liquid 


J  -  A/cm2 


(a) 


water  flooding,  especially  for  the  GDL  beneath  the  bipolar  plate 
ribs. 

The  potential  distribution  inside  the  stack  at  4  A  discharge 
is  shown  in  Fig.  6.  The  distribution  information  is  given  along 


J  -  A/cm2 


(b) 

J  -  A/cm2 


Fig.  9.  Current  density  (A  cm  2)  in  the  membrane:  (a)  cell  1;  (b)  cell  3;  (c)  cell  6. 
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the  intersection  line  of  vz  and  yz  section  planes  across  the  cen¬ 
tral  point  of  the  computational  domain.  Potential  of  the  anode 
side  end  plate  is  set  to  zero  and  as  a  reference  potential,  then 
the  potential  and  overpotential  are  given  according  to  the  ref¬ 
erence  potential.  The  positions  where  the  potentials  jump-up 


correspond  to  the  positions  of  MEAs  in  the  stack.  Activation 
overpotential  in  the  catalyst  layers  (especially  in  cathode  cata¬ 
lyst  layers)  and  ionic  impedance  in  PEM  are  the  main  causes  of 
overpotential  across  the  stack.  The  ohmic  overpotential  in  the 
GDLs  and  bipolar  plates  is  much  smaller  than  the  former  two. 


C  H20  -  mol/mA3 
1.45 

1.45 

1.4 

1.35 

1.3 

1.25  - 

1.2  - 

1.15  — I — 

1.15 


C  02  -  mol/mA3 
8.2 

8.18 
8.16 
8.14 
8.12 
8.1 

8.08  -j 

8.06  - 

8.04  - 

8.02  - 

8 - 

8 


(b) 


Fig.  10.  Hydrogen  (a)  and  water  (b)  molar  concentrations  in  anode  flow  field  of 
cell  3. 


Fig.  11.  Oxygen  (a)  and  water  (b)  molar  concentrations  in  cathode  flow  field  of 
cell  3. 
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With  the  solid  phase  conductivity  data  listed  in  Table  1 ,  it  is  esti¬ 
mated  that  the  electric  resistance  of  the  stack  is  about  30  mQ  and 
ohmic  overpotential  is  120  mV  at  4  A  discharge.  From  the  sim¬ 
ulated  polarization  curve  shown  in  Fig.  5,  the  stack  resistance 
could  be  calculated  with  a  nonlinear  fit  of  the  Tafel  equation 
shown  below: 


ductivity  seems  to  increase  along  x  direction  because  membrane 
conductivity  is  higher  under  higher  temperature  according  to  Eq. 
(9)  [2],  where  o  is  membrane  conductivity  and  X  is  membrane 
water  content: 


o'  =  (0.005 139A 


0.00326)  exp 


E  =  E0  +  Blogi  -  R  i  (8) 

The  overall  stack  resistance  is  fitted  to  be  256  m^,  which 
reveals  that  the  membrane  ionic  resistance  is  about  seven  times 
greater  than  the  electric  resistance.  The  electric  overpotential, 
ionic  overpotential  and  activation  overpotential  are  calculated 
to  be  0. 12  V,  1 .023  V  and  2.285  V,  respectively.  Activation  takes 
most  part  of  the  stack  overpotential. 

Fig.  7(a  and  b)  shows  the  air  velocity  and  total  pressure  distri¬ 
bution  in  the  stack,  respectively.  The  air  inlet  velocity  is  5  m  s_1 
and  the  stack  current  is  4  A.  It  could  be  seen  that  the  air  flow 
rate  increases  along  the  ducts,  which  is  because  the  water  pro¬ 
duction  and  temperature  increase  along  the  flow  direction.  The 
highest  speed  could  be  16.87  m  s-1 .  The  total  pressure  distribu¬ 
tion  is  shown  in  plot  (b).  It  shows  that  the  total  pressure  loss 
is  3 12.4  Pa,  which  includes  294.6  Pa  static  pressure  loss  and 
17.8  Pa  dynamic  pressure  loss  (pu2/2).  The  dynamic  pressure 
loss  is  much  smaller  than  the  static  pressure  loss  because  of  the 
low  density  of  air. 

The  main  issue  for  an  air  cooling  stack  is  not  the  pressure  loss 
but  the  effectiveness  of  stack  cooling.  With  higher  stack  temper¬ 
ature,  water  is  easier  to  evaporate  because  the  saturation  vapor 
pressure  of  water  is  higher.  The  membrane  will  be  inclined  to 
dehydrate  and  lose  ionic  conductivity,  which  will  make  the  fuel 
cell  performance  deteriorate.  In  fact,  air  cooling  stacks  with¬ 
out  external  humidification  could  not  be  steadily  operated  when 
stack  temperature  is  above  323  K  [24].  Fig.  8  shows  temperature 
distribution  in  the  fuel  cell  stack  with  5  m  s_1  air  inlet  velocity 
and  4  A  discharge.  Plot  (a)  shows  temperature  in  yz  and  xz  sec¬ 
tion  planes  and  (b)  in  xy  section  plane  across  the  air  flow  field 
of  cell  3  (numbered  from  the  anode  side).  The  plots  show  small 
temperature  gradient  in  y  direction,  but  great  gradient  in  x  and  z 
directions.  Temperature  in  ME  As  is  higher  than  that  in  the  flow 
field,  which  could  be  seen  clearer  in  Fig.  13.  For  4  A  discharge 
and  5 ms-1  air  inlet  velocity,  air  temperature  increases  about 
5  K  from  the  air  inlet  to  the  air  outlet. 

Current  density  distributions  in  the  membrane  of  cells  1,  3 
and  6  are  shown  in  Fig.  9.  Although  current  density  distributions 
show  some  small  difference  in  different  cells,  the  distribution 
trend  is  similar  in  every  cell:  current  density  decreases  from  the 
air  inlet  to  the  outlet,  and  increases  to  some  extent  from  hydro¬ 
gen  inlet  to  outlet.  The  highest  current  density  regions  appear  on 
the  air  inlet  side  and  away  from  the  hydrogen  inlet  side.  One  may 
doubt  the  distribution  trend  with  some  reasons:  firstly,  oxygen 
reduction  reaction  (ORR)  kinetics  will  increase  along  the  air  flow 
direction  because  cell  temperature  increases  from  air  inlet  to  out¬ 
let  as  shown  in  Fig.  8(b);  secondly,  water  concentration  seems 
to  be  higher  on  the  air  outlet  side  because  of  water  production. 
As  a  result,  it  seems  that  current  density  should  be  higher  on  the 
air  outlet  side  rather  than  the  inlet  side.  Thirdly,  membrane  con¬ 


This  argument  sounds  like  reasonable,  but  there  are  some 
more  important  reasons.  For  the  first  point,  it  should  be  noted  that 
the  influence  of  temperature  to  ORR  kinetics  is  not  considered  in 
the  present  model.  With  Eq.  (10)  as  given  in  [22],  ORR  exchange 
current  density  increases  20%  when  temperature  increases  from 
300  K  to  305  K. 


*0,2  =  *o,i  exp 


A  E  f  1 
~R  \J2 


(10) 


Although  this  omission  will  cause  some  influence,  it  is  not  the 
main  reason  for  the  current  distribution.  Otherwise,  there  would 
not  be  a  great  current  density  gradient  shown  in  the  y  direction 
(Fig.  9),  because  almost  no  temperature  gradient  is  shown  in 
this  direction  (Fig.  8(b)).  The  current  density  gradient  in  the  y 
direction  indicates  that  the  fuel  cell  reaction  is  more  sensitive  to 
humidity  than  ORR  kinetics.  In  Fig.  10,  molar  concentrations  of 
hydrogen  and  water  in  hydrogen  flow  field  of  cell  3  are  shown. 
It  should  be  noted  that  water  concentration  is  highest  in  the  left- 
bottom  corner,  where  the  highest  current  density  is  shown. 

Secondly,  because  the  air  stoichiometry  is  quite  large,  very 
small  oxygen  and  water  concentration  gradients  are  shown.  In 
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Fig.  12.  Membrane  water  content  of  cell  3  at  4  A  discharge. 
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Fig.  13.  Temperature  distribution  of  the  stack  with  different  air  inlet  velocities  in  xz  section  plane:  (a)  5  m  s  1 ;  (b)  3  m  s  1 ;  (c)  1  m  s' 


Fig.  11,  oxygen  and  water  concentrations  in  air  flow  field  of  cell 
3  are  shown.  With  5 ms-1  air  inlet  velocity,  air  flux  is  calcu¬ 
lated  to  be  274  L  min- 1 ,  while  the  stoichiometric  air  flux  for  4  A 
operation  is  434  mL  min-1,  which  indicates  that  the  air  is  630 
times  in  excess.  Although  with  water  production  from  the  elec¬ 
trochemical  reaction,  the  water  concentration  is  actually  shown 
to  decrease  in  the  v  direction  because  of  the  gas  density  decrease 
with  the  increase  of  temperature. 

Thirdly,  it  is  true  that  temperature  will  influence  the  mem¬ 
brane  conductivity,  but  it  will  influence  the  membrane  water 
content  A  greater.  Although  water  concentration  gradients  are 
not  great  in  Figs.  10(b)  and  1 1(b),  water  humidity  will  decrease 
in  v  direction  with  the  increase  of  temperature.  When  temper¬ 
ature  increases  from  300  K  to  305  K,  air  relative  humidity  will 
decrease  from  80%  to  60.2%  because  of  the  32.7%  increase 
of  saturated  vapor  pressure.  So  water  activity  on  both  anode 
and  cathode  sides  will  vary  quite  much  in  v  direction  accord¬ 
ing  to  local  temperatures,  which  will  in  turn  impact  the  mem- 


Fig.  14.  Voltage  and  highest  temperature  in  stack  at  different  air  inlet  velocities 
at  4  A  discharge. 
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brane  water  content.  Fig.  12  shows  the  membrane  water  con¬ 
tent  for  5 ms-1  air  inlet  velocity  and  4 A  discharge.  A  large 
gradient  of  membrane  water  content  is  shown  in  this  plot.  It 
could  also  be  seen  that  the  water  content  distribution  is  very 


z  x 

(b)  - 1 

Fig.  15.  Membrane  water  content  of  cell  3  at  4  A  discharge.  Air  inlet  velocity: 
(a)  2  m  s- 1  and  (b)  1ms-1. 


similar  to  that  of  the  current  density  distribution  shown  in 
Fig.  9(b). 

With  the  premise  of  effective  stack  cooling,  a  lower  air 
flow  rate  will  reduce  the  power  consumption  of  the  blower  and 
increase  the  system  efficiency.  Different  air  flow  rate  cases  were 
studied  to  determine  the  effectiveness  of  cooling.  The  air  inlet 
velocity  is  specified  to  be  5-1  ms-1  and  results  are  shown  in 
Figs.  13  and  14.  In  Fig.  13,  temperature  distributions  in  vz  pro¬ 
file  through  the  central  point  are  shown  in  contour  plots  with 
air  inlet  velocity  5 ms-1,  3 ms-1  and  1ms-1,  respectively. 
From  these  plots  it  could  be  seen  that  the  isothermal  curves 
form  peaks  at  the  positions  of  MEAs  and  large  temperature 
gradients  are  shown  in  the  membranes.  In  the  cathode  catalyst 
layers,  water  is  produced  and  releases  heat,  so  the  tempera¬ 
tures  in  the  cathode  catalyst  layers  are  the  highest  at  the  same 
position  along  the  flow  direction.  The  membranes  have  the  low¬ 
est  thermal  conductivity  and  large  temperature  differences  are 
shown  on  double  sides  of  the  membranes.  The  air  inlet  velocity 
decreases  from  5ms-1  to  3ms-1,  the  highest  stack  temper¬ 
ature  increases  from  306  K  to  310  K.  While  the  inlet  velocity 
decreases  to  1ms-1,  the  highest  temperature  increases  to  about 
330  K.  In  Fig.  14,  the  highest  stack  temperatures  and  stack  volt¬ 
ages  are  shown  according  to  different  air  inlet  velocities.  For 
the  case  of  1ms-1,  the  air  flux  is  too  small  for  the  stack  cool¬ 
ing  and  the  stack  temperature  is  higher  than  323  K.  For  the  case 
of  2ms-1,  the  air  is  about  250  times  in  excess  and  the  stack 
temperature  is  less  than  323  K.  Fig.  15  shows  membrane  water 
content  of  cell  3  with  2ms-1  and  1ms-1  air  inlet  velocity 
at  4  A  discharge.  It  could  be  seen  that  the  water  content  is  in 
a  range  5-9;  but  when  inlet  velocity  decreases  to  be  1ms-1, 
membrane  water  content  decreases  to  very  low,  less  than  4.5. 
At  this  time,  the  fuel  cell  membrane  will  dry  out  and  could  not 
operate  steadily.  In  experiments  the  fuel  cell  stack  performance 
will  decrease  sharply  when  the  air  flow  rate  is  not  enough.  For 
2ms-1  air  inlet  velocity,  although  the  stack  power  is  relative 
low,  11-13. 7  W  at  5ms-1  inlet  velocity  is  the  choice  for  stack 
operation.  It  needs  to  be  pointed  out  that,  because  the  cell  tem¬ 
perature  varies  with  stack  current,  for  a  better  stack  temperature 
management,  it  is  better  to  adjust  the  air  flow  rate  according  to  the 
load. 

5.  Conclusion 

A  mini  PEMFC  stack  was  simulated  three  dimensionally  with 
a  reaction  and  transport  coupled  CFCD  model.  With  reasonable 
simplification,  gridpoints  of  the  computational  domain  could  be 
effectively  reduced  and  allow  a  full  3D  fuel  cell  stack  compu¬ 
tation  to  be  conducted  in  personal  computers.  The  simulation 
results  were  compared  with  experimental  results  and  reasonable 
agreement  was  found.  Temperature  distributions  in  the  stack 
with  different  air  flow  rates  were  given  and  the  results  show 
that,  for  4  A  operation,  the  air  inlet  velocity  could  be  as  low 
as  2ms-1  (250  times  excess)  and  allow  the  stack  to  be  effec¬ 
tively  cooled,  while  for  1ms-1,  the  membrane  water  content 
was  too  low  to  keep  the  stack  in  steady  operation.  The  simu¬ 
lation  results  illustrate  a  method  for  large-scale  fuel  cell  stack 
simulation. 
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